package fem2;

import inf.jlinalg.Array2DMatrix;
import math2.ConstantFunctionRnToR;
import math2.FunctionRnToR;
import math2.PulseFunctionRToR;

import org.junit.Test;

import fem2.analysis.Analysis;
import fem2.analysis.FreeVibrationAnalysis;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.MassMatrixType;
import fem2.enu.SolverType;
import fem2.enu.State;
import fem2.material.StVenantKirchhoffMM;
import fem2.observer.TimeStepObserver;
import fem2.strategies.Strategy;

public class PatchTestDynamics extends AbstractStructuralStaticDemo {

	private int option = 1;

	public PatchTestDynamics() {
		projectDir = "/home/hbui/kratos_janosch";
		projectName = "patchtest-dynamics.gid";
		meshFileName = "patch-freevib";
	}

	@Override
	public Mesh createMesh() {
		if (option == 1)
			return createMeshQ4();
		else if (option == 0)
			return createMeshT3();
		else if (option == 2)
			return createMeshT3_2();
		else if (option == 3)
			return createMeshQ4_2();
		else
			throw new InvalidOptionError();
	}

	public Mesh createMeshT3() {
		Mesh mesh = new Mesh();

		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(1, 0));
		Node n3 = mesh.addNode(new Node(1, 1));
		Node n4 = mesh.addNode(new Node(0, 1));

		mesh.addFace(new Face(n1, n2, n3));
		mesh.addFace(new Face(n1, n3, n4));

		mesh.setNumberOfFields(2);

		n1.setConstraint(new Constraint(false, false));
		n4.setConstraint(new Constraint(false, true));
		n2.setConstraint(new Constraint(true, false));

		return mesh;
	}

	public Mesh createMeshT3_2() {
		Mesh mesh = new Mesh();

		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(1, 0));
		Node n3 = mesh.addNode(new Node(1, 1));
		Node n4 = mesh.addNode(new Node(0, 1));
		Node n5 = mesh.addNode(new Node(0.5, 0.5));

		mesh.addFace(new Face(n1, n2, n5));
		mesh.addFace(new Face(n5, n2, n3));
		mesh.addFace(new Face(n5, n3, n4));
		mesh.addFace(new Face(n5, n4, n1));

		mesh.setNumberOfFields(2);

		n1.setConstraint(new Constraint(false, false));
		n4.setConstraint(new Constraint(false, true));
		n2.setConstraint(new Constraint(true, false));

		return mesh;
	}

	public Mesh createMeshQ4() {
		Mesh mesh = new Mesh();

		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(1, 0));
		Node n3 = mesh.addNode(new Node(1, 1));
		Node n4 = mesh.addNode(new Node(0, 1));

		mesh.addFace(new Face(n1, n2, n3, n4));

		mesh.setNumberOfFields(2);

		n1.setConstraint(new Constraint(false, false));
		n4.setConstraint(new Constraint(false, true));
		n2.setConstraint(new Constraint(true, false));

		return mesh;
	}

	public Mesh createMeshQ4_2() {
		Mesh mesh = new Mesh();

		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(1, 0));
		Node n3 = mesh.addNode(new Node(1, 1));
		Node n4 = mesh.addNode(new Node(0, 1));
		Node n5 = mesh.addNode(new Node(0.5, 0));
		Node n6 = mesh.addNode(new Node(0.5, 1));

		mesh.addFace(new Face(n1, n5, n6, n4));
		mesh.addFace(new Face(n5, n2, n3, n6));

		mesh.setNumberOfFields(2);

		n1.setConstraint(new Constraint(false, false));
		n4.setConstraint(new Constraint(false, true));
		n2.setConstraint(new Constraint(true, false));
		n5.setConstraint(new Constraint(true, false));

		return mesh;
	}

	@Override
	public MaterialModel createMaterial(Model m) {
		double E = 2;
		double nu = 0.3;
		double thickness = 1.0;
		double rho = 1.0;
		State ss = State.PLANE_STRESS;

		return new StVenantKirchhoffMM(E, nu, thickness, rho, ss);
	}

	@Override
	public Element createElement(MeshPart mp, MaterialModel mm) {
		return new StructuralElement(mp, mm);
	}

	@Override
	public Model createConditions(Model m) {
		Node n1 = MeshUtilities.seekNode(m.getMesh(), 1, 0);
		Node n2 = MeshUtilities.seekNode(m.getMesh(), 1, 1);
		Edge e = m.getMesh().addEdge(new Edge(n1, n2));
		e.setNumberOfFields(2);

		FunctionRnToR lambda = new PulseFunctionRToR(0, 1e-6, 1);
		FunctionRnToR p = new ConstantFunctionRnToR(1.0);

		Load f = new DistributedLoad(e, p, null);
		f.setTimeFactor(lambda);
		m.addLoad(f);

		return m;
	}

	@Test
	public void test1T3() {
		option = 0;
		runDynamics(MassMatrixType.DIRECT_MASS_LUMPING);
	}

	@Test
	public void test2T3() {
		option = 0;
		runDynamics(MassMatrixType.VARIATIONAL_MASS_LUMPING);
	}

	@Test
	public void test2T3_2() {
		option = 2;
		runDynamics(MassMatrixType.VARIATIONAL_MASS_LUMPING);
	}

	@Test
	public void test1Q4() {
		option = 1;
		runDynamics(MassMatrixType.DIRECT_MASS_LUMPING);
	}

	@Test
	public void test2Q4() {
		option = 1;
		runDynamics(MassMatrixType.VARIATIONAL_MASS_LUMPING);
	}

	@Test
	public void test1Q4_2() {
		option = 3;
		runDynamics(MassMatrixType.VARIATIONAL_MASS_LUMPING);
	}

	/**
	 * 
	 */
	private void runDynamics(MassMatrixType massType) {
		Model m = createModel();

		for (int i = 0; i < m.countElements(); i++) {
			int n = m.getElement(i).getMeshPart().countDofs();
			Array2DMatrix me = new Array2DMatrix(n, n);
			m.getElement(i).computeMe(me, massType);
			// m.getElement(i).computeMe(me,
			// MassMatrixType.DIRECT_MASS_LUMPING);
			Debugger.watch(me);
		}

		Analysis an = new FreeVibrationAnalysis(m, massType, 0.0, 0.0); // no
																		// damping
																		// ((FreeVibrationAnalysis)
																		// an).setAnalysisTypeToNonAdaptiveNewmark(NewmarkType.NewmarkAlpha,
		// 0.8, 0.0, 0.1, 10);
		((FreeVibrationAnalysis) an).setAnalysisTypeToCDM(0.0, 0.1, 10);

		Strategy s = an.getStrategy();
		// an.getStrategy().setEchoLevel(EchoLevelType.AUX2,
		// EchoLevelType.AUX3);
		s.setEchoLevel(EchoLevelType.OUTPUT, EchoLevelType.DEBUG);
		s.setSolverType(SolverType.BAITSCH_SOLVER_GENERIC);

		// Node n2 = m.getMesh().getNode(1);
		// Node n3 = m.getMesh().getNode(2);
		// Node n4 = m.getMesh().getNode(3);
		// s.addObserver(new DisplacementObserver(n2, n3, n4));

		an.run();
	}

	@Override
	public void run() {
		option = 3;
		Model m = createModel();

		Analysis an = new FreeVibrationAnalysis(m, MassMatrixType.VARIATIONAL_MASS_LUMPING, 0.0,
				0.0); // no damping
		// ((FreeVibrationAnalysis)
		// an).setAnalysisTypeToNonAdaptiveNewmark(NewmarkType.NewmarkAlpha,
		// 0.8, 0.0, 0.1, 10);
		((FreeVibrationAnalysis) an).setAnalysisTypeToCDM(0.0, 0.1, 10);

		Strategy s = an.getStrategy();
		s.setSolverType(SolverType.BAITSCH_SOLVER_GENERIC);
		s.setEchoLevel(EchoLevelType.OUTPUT, EchoLevelType.AUX1);

		addObservers(m, s, new TimeStepObserver(s));

		an.run();

		// postProcess(m);

		Debugger.warn("Analysis completed");
	}

	public static void main(String[] args) {
		AbstractStructuralStaticDemo demo = new PatchTestDynamics();

		demo.run();
	}
}
